library(dplyr)
library(ggplot2)
library(visreg)
library(texreg)
library(lme4)
library(stargazer)
library(interactions)

setwd("~/Location/of/Data")
# Load member-level data
members <- read.csv("member_level_data.csv")

# load trip-level data
trips <- read.csv("trip_level_data.csv")

#### DESCRIPTIVES ####

# Calculating total number of days traveled
sum(members$days.traveled, na.rm = T) # total days traveled by all members over study period
sum(members$days.traveled)/sum(!is.na(members$days.traveled)) # average number of days per trip

# Table 1
# # Tables of number of trips per member for different types of member
# committee membership, chair status, seniority, extremity, partisanship
party <- members %>% group_by(democrat) %>%
   summarize(trip.avg = mean(trips),
             trip.sd = sd(trips),
             trip.min = min(trips),
             trip.max = max(trips),
             n = n())

chair <- members %>% group_by(chair) %>%
   summarize(trip.avg = mean(trips),
             trip.sd = sd(trips),
             trip.min = min(trips),
             trip.max = max(trips),
             n = n())

power <- members %>% group_by(power.committee) %>%
   summarize(trip.avg = mean(trips),
             trip.sd = sd(trips),
             trip.min = min(trips),
             trip.max = max(trips),
             n = n())

members$seniority3 <- cut_number(members$seniority, 3, ordered_result = T)
senior <- members %>% group_by(seniority3) %>%
   summarize(trip.avg = mean(trips),
             trip.sd = sd(trips),
             trip.min = min(trips),
             trip.max = max(trips),
             n = n())

members$extremity3 <- cut_number(members$extremity, 3, ordered_result = T)
extremity <- members %>% group_by(extremity3) %>%
   summarize(trip.avg = mean(trips),
             trip.sd = sd(trips, na.rm = T),
             trip.min = min(trips),
             trip.max = max(trips),
             n = n()) %>%
   filter(!is.na(extremity3))

leader <- members %>% group_by(leader) %>%
   summarize(trip.avg = mean(trips),
             trip.sd = sd(trips),
             trip.min = min(trips),
             trip.max = max(trips),
             n = n())

members$winpct3 <- cut_number(members$winpct, 3, ordered_result = T)
winpct <- members %>% group_by(winpct3) %>%
   summarize(trip.avg = mean(trips),
             trip.sd = sd(trips),
             trip.min = min(trips),
             trip.max = max(trips),
             n = n()) %>%
   filter(!is.na(winpct3))

des <- bind_rows(party, leader, chair, power, senior, winpct, extremity)

Variable <- c("Republican", "Democrat", 
              "Rank and File", "Leader",
              "Non-Chair", "Chair",
              "Non-Power Committee", "Power Committee",
              "Low Seniority", "Mid Seniority", "High Seniority",
              "Low Win Pct.", "Mid Win Pct.", "High Win Pct.",
              "Low Extremity", "Mid Extremity", "High Extremity")

Mean <- des$trip.avg
SD <- des$trip.sd
Mins <- des$trip.min
Maxes <- des$trip.max
N <- des$n

# table 1 output
stargazer::stargazer(data.frame(Variable, Mean, SD, N), 
                     digits = 1, summary = F)

# Histogram of trips per congress
# Figure 1
ggplot(members, aes(trips)) + geom_histogram(binwidth = 1) + theme_classic() +
      xlab("Trips") +
      ylab("Number of Members")
#ggsave("figure1.png", width = 6, height = 6*.66)

## Total trips per congress over time
# Figure 2
trips %>% group_by(cong) %>% summarize(tot = n()) %>%
      ggplot(aes(cong, tot)) + geom_col(aes(group = factor(cong))) + 
      labs(x="Congress", y="Total Trips per Congress") +
     # ggtitle("Total Gift Trips per Congress") +
      theme_classic() + guides(fill="none")


# count of trips by member over time
# Figure 3
trips %>% ggplot(aes(factor(MemberTrip))) + 
      geom_bar(aes(fill = factor(cong)), position = "dodge") + 
      labs(x="Congresses (110th to 115th)", y="Total Trips per Congress") + 
      #ggtitle("Total Gift Trips per Congress by Staffers and Members") +
      scale_x_discrete(labels=c("0" = "Staffers", "1" = "Members")) + 
      scale_fill_grey() + theme_classic() + guides(fill = "none")


# % of group trips that are overseas
trips %>% filter(n_travelers > 1) %>%
      group_by(TripID, Foreign) %>%
      summarize(n = n()) %>%
      ungroup() %>%
      summarize(perc.foreign = mean(Foreign))

# % of solo trips that are overseas
trips %>% filter(n_travelers == 1) %>%
      group_by(TripID, Foreign) %>%
      summarize(n = n()) %>%
      ungroup() %>%
      summarize(perc.foreign = mean(Foreign))


#### DV: TRIPS, NEGATIVE BINOMIAL MODEL ####


### RANDOM-EFFECTS MODELS

# Table 2, Model 1
trips.re1 <- glmer.nb(trips ~ winpct + chair + subchair + power.committee + leader + seniority + 
                         extremity * democrat + female + black + latino + LES +
                         (1|ICPSRnumber) + (1|cong), data = members)
screenreg(list(trips.re1), digits = 3)


### PREDICTED TRIPS ###
# party leader effect
visleader <- visreg(trips.re1, xvar = "leader", scale = "response", rug = F,  overlay = T)
range(visleader$fit$visregFit)[2] - range(visleader$fit$visregFit)[1]

# committee chair effect
vischair <- visreg(trips.re1, xvar = "chair", scale = "response", rug = F,  overlay = T)
range(vischair$fit$visregFit)[2] - range(vischair$fit$visregFit)[1]

# win pct effects
viswinpct1 <- visreg(trips.re1, xvar = "winpct", scale = "response", rug = F,  overlay = T)
comp1 <- filter(viswinpct1$fit, winpct >= 51 & winpct <= 80)
range(comp1$visregFit)[2] - range(comp1$visregFit)[1]

comp2 <- filter(viswinpct1$fit, winpct >= 51 & winpct < 101)
range(comp2$visregFit)[2] - range(comp2$visregFit)[1]

# ideological extremity effect
# for Dems
visext <- visreg(trips.re1, xvar = "extremity", by = "democrat", 
                 scale = "response", rug = F,  overlay = T)
demext <- visext$fit %>% 
   filter(democrat == 1)
range(demext$visregFit)[2] - range(demext$visregFit)[1]

# for Republicans
repext <- visext$fit %>% 
   filter(democrat == 0)
range(repext$visregFit)[2] - range(repext$visregFit)[1]


#covariates set to non-leader values, diff = 2.39
visleader <- visreg(trips.re1, xvar = "leader", scale = "response", rug = F,  overlay = T,
                    cond = list(seniority = 5.566, extremity = .459))
#covs set to leader values, diff = 2.37
visleader <- visreg(trips.re1, xvar = "leader", scale = "response", rug = F,  overlay = T,
                    cond = list(seniority = 7.202, extremity = .542))

# covs set to non-chair values, diff = 4.09
vischair <- visreg(trips.re1, xvar = "chair", scale = "response", rug = F,  overlay = T,
                   cond = list(seniority = 5.416, extremity = .475))
# covs set to chair values, diff = 4.43
vischair <- visreg(trips.re1, xvar = "chair", scale = "response", rug = F,  overlay = T,
                   cond = list(seniority = 10.452, extremity = .237))


#### DV: LES, LINEAR MODEL ESTIMATION ####

## Random-Effects Models 

# Main models
# Table 2, Model 2
les.re1 <- lmer(LES ~ trips +  winpct + chair + subchair + power.committee + leader + seniority + 
                      extremity + democrat + female + black + latino + statelegxprof +
                      (1|ICPSRnumber) + (1|cong), data = members)

lesvis <- visreg(les.re1, xvar = "trips", partial = F, xlim = c(0,20))
range(lesvis$fit$visregFit)[2] - range(lesvis$fit$visregFit)[1]


# Table 2, Model 3
les.re2 <- lmer(LES ~ trips*democrat + winpct + chair + subchair + power.committee + leader + seniority + 
                      extremity + democrat + female + black + latino + statelegxprof +
                      (1|ICPSRnumber) + (1|cong), data = members)
# Table 2
screenreg(l = list(les.re1, les.re2), digits = 2)

# Figure 4
lesvis2 <- interact_plot(les.re2, trips, democrat, colors = c("black", "black"),
                         modx.labels = c("Republican", "Democrat"), interval = T,
                         legend.main = "Party", x.label = "Trips")
lesvis2



#### Special-trip models 

# Solo vs. group trips

les.re9 <- lmer(LES ~ trips_solo + trips_group + winpct + chair + subchair + power.committee + leader + seniority + 
                      extremity + democrat + female + black + latino + statelegxprof +
                      (1|ICPSRnumber) + (1|cong), data = members)

# Table 4
screenreg(l = list(les.re9), single.row = T, digits = 3)


# Domestic vs. foreign
# Table 6
les.re3 <- lmer(LES ~ domestic.trips + foreign.trips + winpct + chair + subchair + power.committee + leader + seniority + 
                      extremity + democrat + female + black + latino + statelegxprof +
                      (1|ICPSRnumber) + (1|cong), data = members)

les.re4 <- lmer(LES ~ domestic.trips*democrat + foreign.trips*democrat + winpct + chair + subchair + power.committee + leader + seniority + 
                      extremity + democrat + female + black + latino + statelegxprof +
                      (1|ICPSRnumber) + (1|cong), data = members)
# Table 6
screenreg(l = list(les.re3, les.re4), digits = 2)



# Staffer vs. members
les.re5 <- lmer(LES ~ staffer.trips + member.trips + winpct + chair + subchair + power.committee + leader + seniority + 
                      extremity + democrat + female + black + latino + statelegxprof +
                      (1|ICPSRnumber) + (1|cong), data = members)
# Table 7
screenreg(les.re5, single.row = T, digits = 3)

# average number of guests per event
gpe <- trips %>% group_by(TripID) %>%
   summarize(guestsPerEvent = n())
summary(gpe$guestsPerEvent)
proportions(table(gpe$guestsPerEvent))


##### SPONSOR-LEVEL ANALYSES #####
sponsors <- read.csv("travel_sponsor_level_data.csv")

# Table 5
# Rank sponsors by total.guests hosted in data set
top <- sponsors %>% group_by(TravelSponsor) %>%
   summarize(total.guests = sum(guests),
             total.events = sum(events),
             avg.guestsPerEvent = mean(guestsPerEvent)) %>%
   mutate(guest.prop = total.guests / sum(total.guests)) %>%
   arrange(desc(total.guests))
top$cum.prop <- cumsum(top$guest.prop)

top <- left_join(top, sponsors[, c("TravelSponsor", "sponsor.type")],
                 by = "TravelSponsor")
top <- distinct(top)
top25 <- top[1:25,]
top25 <- mutate(top25, guest.prop = round(guest.prop, 2),
                cum.prop = round(cum.prop, 2),
                avg.guestsPerEvent = round(avg.guestsPerEvent, 0))
top25 

# Count of guests by sponsor type
sponsors %>% group_by(sponsor.type) %>%
   summarize(guestsByType = (sum(guests, na.rm=T))) %>%
   mutate(propOfGuests = guestsByType / sum(guestsByType, na.rm = T)) %>%
   ggplot(aes(sponsor.type, propOfGuests)) + geom_col()


## TABLE 8
# DV = guests
guests.re2 <- glmer.nb(guests ~ lobbyist.hired +
                          (1|TravelSponsor) + (1|cong), data = sponsors)

# DV = events
events.re2 <- glmer.nb(events ~ lobbyist.hired +
                          (1|TravelSponsor) + (1|cong), data = sponsors)
screenreg(list(events.re2, guests.re2), single.row = T, digits = 3)

